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In this work, we present a new fitting of the Na + HF potential energy surface (PES) utilizing 
a new optmization method based in Genetic Algorithm. Topology studies, such as isoenergetic 
contours and Minimum Energy Path (MEP), show that the quality of this new PES is comparable 
to the best PES of literature. These facts, suggests that this new approach can be utilized as new 
tool to fit PES of reactive systems. 



I. INTRODUCTION 

Genetic algorithms Q, |^ |^ ^ IE 1^ have been applied successfully in the description of a variety of global mini- 
mization problems. It have as well attracted significant attention due to their suitability for large-scale optimization 
problems, specially for those in which a desired global minimum is hidden among many local minima. 

The main object of this paper is to propose a genetic algorithm optimization technique(GAOT) for fitting the 
PES through electronic energies obtained by ab initio calculations. In order to present and to test the method, we 
reproduce the PES of the reaction involving HE and an alkali metal, namely 

Na{?,'^Si/2) + HF{X^T.+) ^ NaF{X^j:+) + HCSi/2), (1) 

utilizing as a trial function a Bond Order (BO) polynomial expansion 01 as well as ab initio calculations published by 
Lag ana et al^ . The motivation behind the choice of this reaction is its high endoergicity and its bent transition state 
|l3l |. Further more, this reactive process has been experimentally theoretically 0,0] and computational 

investigated [T^ . 

This paper is organized as follows. In Section 2, we present the main characteristics of the GAOT. The details of 
the GAOT fitting and its comparison with other SEPs are shown in the Section 3. Our conclusions are contained in 
the Section 4. 

II. MODEL 
A. The Problem 

In order to fit a given functional form y([a],r) in some set of Up points (rp,ep), we want that the GAOT finds a 
set of parameters [a] = [ai, a2, .., am] that minimize the mean square deviation 

Tip Tip 

s = Y.^l^Y.i^p--^pf (2) 

p p 

where Cp = V'([a],rp). 

B. Codification 



In our genetic algorithm the population is coded in a binary discrete cube named A, with I x m x n bits. The 
elements of A, a^Tc, are either or 1, with i, j, k integers numbers 1 < i < n, I < j < m, 1 < k < n. The label i refers 
to the component i of the gene j of the individual k. Therefore, A represents a population of n individuals, each one 
of them have a genetic code with m genes. Each gene is a binary string with I bits. 



* roncaratti@fis.unb.br 



2 



The genetic code of the individual k is given by 

[a]/j = \aik,a2k,--,amk], 

were 

^jk = ^2*~^ayfe (3) 

i=l 

is a integer number composed with the binary string aijka2jk--0-ijk--0'i3k- It is defined on the interval [0,2' — 1]. To 
define the real search space for each parameter, we transform 

ajk ajk = ^J~[ + ^3 

were Ujk is a real number defined on the interval dj = [a™*", a™""^]. 

<5,- = [a7™,a7'*^]. (5) 
Now we define the phenotype of the individual k,Vk = F([a]fc, r) where 

[a]fe = [aik, a2k, ajk, ■ ■, amk] (6) 
is a set of coefiicients that characterize the individual k. With this we define the fitness of a phenotype k 

where 

Sk = J2(^kpf = J2(''p-Vk^f (7) 
p p 

and Smax is worst individual in the population. Vkp = V{[a]k,i%) and Skp is the difference among the ab initio energy 
Bp and the fit of the individual k in the configuration f = fp. 

C. Operators 

We use the most common operators: selection, recombination and mutation. The selection operator normalize the 

vector Sk 

that represents the probability of each individual been selected for a recombination through a roulette spinning. For 
the purpose of this work we selected n/2 individuals (parents) that will generate, through the recombination operator, 
n/2 new individuals (offsprings). So, to make a new generation we joint the n/2 old strings (parents) with a, n/2 
new strings (offsprings) in order to maintain a population with fixed number n. The recombination operator is a 
cross-over operator that recombine the binary string of each gene j of two random selected individuals to form two 
new individuals. In this work we use a two random point cross-over. 

The mutation operator flip Nmut random selected bits in a population. We choose Nmut to make the probability 
of change of a given bit equal to 0.01 per cent. So, in a population of I x mx n bits, we make 



I xmx n 

where q is the probability of change in one bit. 

The elitist strategy consists of copying an arbitrary number N^i of the best individual on the population in the 
next generation. It warrants that this individual will not be extinguished. 
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III. FITTING THE AB INITIO PES 



Lagana et al. construed the BOS PES of the Na+HF reaction considering a total of the 425 ah initio values 8] , being 
42 values (Table 2 of the Ref.Q) calculated in the region that better characterize the collinear Na-HF and F-NaH 
geometries (insertion). All these 425 energy values cover a relevant portion of the surface at 9, the angle formed by 
the NaF and HF internuclear distances, equal to 6'=0°, 45°, 60°, 75°, 90°, 120° and 180°. Analytical representations 
of the B05 PES were obtained using a BO polynomial expansion for both two- and three-body termsQ, following the 
standard many-body form. Each two-body term was construed fitting a polynomial of the fourth order in the related 
BO variables. The three-body term was fitted using the BO polynomial expansion given by 

5 5 5 

V{RNaF, RhF, RNan) = ^ ^ ^ ^^xyzT]^ aF^H f'^In aH 

(10) 

with X + y + z < 6 and at least two indices differing from zero. The quantities 77™ are defined as 

mi - e-P"^^'"-^-'''^ (11) 

with M=NaF, HF and NaH, where the parameters values of the Na-fHF reactive process are ReNaF — 1.92595 A, 
ReHF = 0.91681 A, ReNaH = 1.88740 A, PNaF = 0.88260 A'\ Phf = 2.19406 A"! and ^NaH = 1.19798 A'^. 

In order to use the GAOT method to reproduce the Na-|-HF PES, we used only the BO polynomial expansion for 
three-body term given by Eq. H10|) . with the same powers utilized to fit the B05 PES. It should be pointed that we 
used only 243 of a total of the 425 ab initio values used to produce the B05 PES. This number was the same used in 
the GSA PES fitting [13 (see table table 1 of the Ref.|3). 

Using the notation of Section II and Eq. (|10l) we define the phenotype of the individual k 



m 

Vk = V{[a]k, RngF, RhF, RNan) = ^ O'jkVNaFVnF ^^O-NaH (12) 

i=i 

where m is the number of coefficients that we use in the expansion. Each coefficient Uj^ had a specific fixed 
combination of powers {xj, ijj, Zj). In this way, we guarantee that the fittest individual in the population had a set of 
coeflficients [a]fc (01 that better fit the Eq. ifTUIl . 

In this work, we define as an acceptable solution the set of 77 coefficients Q that fit the expansion (|10() to 243 
ab initio points in a way that the root mean square deviation |(2Jl be less than 1,0 Kcal/mol. In fact, we can find a 
large number of acceptable solutions. The set of all solutions is the definition of search space (F). The length of F is 
defined by the number m of coefficients and the length I of the binary codification. Each one of the 77 coefficients, 
that define the individual fc, can assume 2' distinct values. So, an individual in the population is only one possibility 
among 2'^™. This value defines the length of F. The precision of F describes the number of digits that are used to 
express a real value ajk and shows the minimal difference between two possible values of Ujk- Being each coefficient 
defined on an arbitrary interval Sj (0), the precision of the coefficient ajk is 



2' 



(13) 



If we do not have any information about the order of magnitude of the Ojk values, we must choose the intervals 
(5j's such that they cover the greatest number of values. However, after some generations, we obtained more precise 
information about the order of magnitude of each coefficient ajk- In order to improve the performance of the standard 
GAOT, we include in our technique the concept of dynamic search space. It consist in the use of information of 
past generations to determine the length and precision of the search space for the next generations. For the first 
generations, when we have few information about the order of magnitude of the coefficients, we do not need many 
digits to represent a real number a^fc, that is, we use a low precision codification given by a low value of /. In this 
way, we make F a "small" search space and the GAOT can find the regions of acceptable solutions faster. Once found 
some of these regions we can redefine the intervals 5j 's and rise the precision rising the length of binary codification 
I. After extensive trials of the parameters values we take m = 77, n — 100, g = 0,01 and N^i = 10. Beside that, 
we always start the GAOT with a random population defined in the initial intervals 5j = [a™™,a™''^] = [—10'"', lO'"'] 
and set the initial value for the length of the binary codification I = 12. In this way we had a search space of length 
2ixm _ 212x77 _ 2^24 q^^^ ^jjg minimal difference of two possible values of a^k is ^ 49. After 1000 generations we 
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redefine I = I + A and = [aj^™, a™"^] where af'" = Qjbest + o™™ x lO"!, aj'"^ = ajbest + x lO"! and ajbest is 
the fitest individual in the population found along the last 1000 generations. We set 10000 generations for each run 
of the GAOT. It should be pointed out that the algorithm is very robust and works properly with an wide range of 
these parameters. 

Coefficients and powers of the polynomial given by Eq^Jfor the GAOT PES are showed in the table |l| We plot in 
the figureClthe GAOT PES, considering 6'=30.0°(a), 6'=180.G°(b) and 61=77.2° (c), with the respectives isoenergetic 
contours. The dashed contours are taken between -160 and -40 kcal/mol and them are spaced each other by 5 kcal/mol. 
One can see that the global shape of the GAOT PES is closed similar to both GSA and B05 PES 17]. 

To better test the new PES, we plot in figure|2]the GAOT fixed angle minimum energy paths (GAOT MEP) of the 
Na-I-HF reaction considering 61=30. 0° (a), 61=180.0° (b) and 9=77.2° (c), respectively. In these figures are also shown 
both GSA and BOS MEPs. In aU these MEPs the zero energy was set at the B05 Na-|-HF asymptote. AU these 
figures show that overall shape of the GAOT, B05 and GSA MEP are very similar. At 6'=77.20° the barrier of the 
B05 reaction is minimum and increases when moving towards collinear or towards more bent geometries. The same 
value was found for GAOT PES. In the table ^ are represented the values of the reactant energy, product energy, 
barrier height and well depth of the B05, GSA and GAOT MEPs considering the following values of the 6'=30°, 60°, 
77.20°, 90°, 120°, 150° and 180°. In the region of the Na-fHF reactant, the differences of the energies find between 
the B05-GA0T and GSA-GAOT MEP, for aU values of the 6'(see theHH), are about 0.02 kcal/mol and 0.73 kcal/mol, 
respectively. However, in the NaF-|-H product region these differences are about 0.14 and 0.48 kcal/mol, respectively. 
In the well region, the maximum differences of the energies between B05-GA0T and GSA-GAOT are about 0.26 
and 0.52 kcal/mol, respectively. The maximum values of these differences in the barrier region are about 0.5 and 1.2 
kcal/mol, respectively. 

IV. CONCLUSIONS 

In this work, we have presented the GAOT method as a new option to fit PES for reactive system. Plots of the 
GAOT PES and B05, GSA and GAOT MEP of the Na-hHF system were made considering different values of the 
NaFH angle. From comparison among these plots, we concluded that these PES have the same global shape. The 
values of the reactant energy, product energy, barrier height and well depth of these MEP at the principal NaFH 
angle considered were very small, within of the error acceptable to reactive process. This ample topologies studies 
reveals that quality of the GAOT PES is comparable the with the best PES find for Na+HF reaction, i.e, B05 and 
GSA PES. 

In a future work, we will present a complete study of the dynamics properties of the Na-|-HF reaction utilizing the 
GAOT PES. These properties will be compared with the dynamics properties determined with the both B05 GSA 
PES. 
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MEP 


Reactant 


Product 


Barrier 


Well 


o 

30 


B05 
GSA 
GAOT 


-0.2303 
-0.9266 
-0.2144 


17.6297 
17.9334 
17.4871 


77.6465 
78.7675 
77.5946 


-1.7699 
-2.0017 
-1.7776 




60 


B05 
GSA 
GAOT 


-0.1978 
-0.9001 
-0.1819 


17.5697 
17.8983 
17.4256 


23.2639 
23.3119 
23.0681 


-3.1729 

-3.0291 
-3.3211 


77.20° 


BOS 

GSA 
GAOT 


-0.1817 

-0.8869 
-0.16S9 


17.5518 
17.8874 
17.4072 


18.6754 

18.5825 
18.5229 


-6.1842 

-5.9927 
-6.3219 
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FIG. 1: The figure shows the isoenergetic contours of the GAOT PES, at 6'=30.0°(a), 6'=180.0°(b) and 6'=77.2°(c). The zero 
energy was set at the Na+HF asymptote. The dashed contours are taken between -160 and -40 kcal/mol and them are spaced 
each other by 5 kcal/mol. 
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FIG. 2: The figure shows the plots of the BOS, GSA and GAOT MEP for the Na+HF reaction calculated at 6 
6 =180.0° (b) and =77.2° (c). The $ angle is associated with the definition of the MEP (see text for discussion). 



=30.0° (a), 



